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Abstract 

The  goal  of  this  work  is  to  achieve  prototype  simulations  of  the  MLEK  process  for 
growth  of  InP  crystals.  We  have  made  several  key  advances  in  this  work,  including 
numerical  simulations  of  three-dimensional  time- dependent  low  Prandt  1-number  melt 
convection,  with  differential  rotation,  and  magnetic  fields.  To  do  this  we  have  developed 
a  3D  general-geometry  spectral  element  code.  We  have  demonstrated  how  rotation  and, 
especially,  differential  rotation  can  stabilize  melt  flows  and  how  magnetic  fields  can  also 
stabilize  these  flows.  We  have  also  shown  that  magnetic  fields,  in  certain  parameter 
regimes,  can  act,  somewhat  surprisingly,  to  destabilize  the  flow.  We  have  also  extended 
our  renormalization-group-based  turbulence  models  to  describe  stratified  flows  across 
the  full  range  of  low-moderate-high  stratification. 
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1  Introduction 


Indium  Phosphide  (InP)  offers  much  promise  for  advanced  photoelectronic  and  electronic 
device  development  because  of  its  attractive  physical  properties  that  allow  high  power,  high 
frequency,  and  short  gate  lengths.  However,  the  growth  of  InP  crystals  involves  a  number  of 
major  technological  challenges  because  of  the  need  to  use  high-pressure  environments  (due 
to  the  high  vapor  pressure  of  the  melt),  the  difficulty  of  working  with  the  volatile  component 
phosphorus,  the  complex  stoichiometry  of  the  mixtures  involved,  and  the  sensitivity  of  the 
grown  crystals  to  defects,  twinning,  and  other  imperfections.  An  attractive  method  to  grow 
InP  crystals  in  practice  is  the  use  of  the  magnetic  liquid  encapsulated  Kyropoulos  (MLEK) 
technique  but  the  complexity  of  this  process  limits  the  opportunity  to  use  trial-and-error 
experimental  techniques  to  optimize  the  crystal  growing  process. 

In  our  work  under  AFOSR  Contract  F49620-94-C-0035,  we  have  made  much  progress 
in  the  simulation  of  the  flow  fields  involved  in  the  MLEK  and  related  growth  processes. 
We  have  made  several  key  advances  in  this  work,  including  numerical  simulations  at  low 
Prandtl  number  melt  convection  with  rotation,  differential  rotation,  and  magnetic  fields.  We 
have  demonstrated  how  rotation  and,  especially,  differential  rotation  can  stabilize  melt  flows 
and  how  magnetic  fields  can  also  stabilize  these  flows.  We  have  also  shown  that  magnetic 
fields,  in  certain  parameter  regimes,  caji  act,  somewhat  surprisingly,  to  destabilize  the  flow. 
We  have  developed  a  simplified  model  that  illustrates  this  effect  of  magnetic  fields.  We 
have  also  extended  our  renormalization-group-based  turbulence  model  to  describe  stratified 
flows  across  the  full  range  of  low-moderate-high  stratification.  Our  simulations  give  the 
first,  spectrally  accurate,  simulations  of  3D  low  Prandtl  number  convection  with  coupled 
rotation,  differential  rotation,  and  magnetic  field.  As  such,  we  believe  that  they  represent  a 
change  in  the  state-of-the-art  of  crystal  growth  modeling  and  simulation.  In  addition,  our 
turbulence  transport  model  for  stratification  should  find  wide  application  to  problems  of 
Air  Force  interest.  For  example,  the  new  model  can  enable  advances  in  the  problem  of  the 
dynamics  of  aircraft  wake  vortices  in  stratified  atmospheres,  a  problem  with  applications  to 
both  civilian  and  military  aircraft  (e.g.,  if  stratification  can  cause  trailing  vortices  to  bounce 
off  a  stratified  atmospheric  layer  into  the  path  of  trailing  aircraft). 
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2  Technical  Background  -  Simulation  of  Melt  Flows 


We  have  studied  flows  in  crystal  melts  in  several  situations.  First,  we  consider  a  model 
problem  consisting  of  a  cylindrical  crucible  of  radius  R  and  height  H  (see  Figure  1).  Here  we 
have  a  prototype  crystal  with  radius  r,  such  that  r/R  =  0.5,  being  pulled  from  the  top  surface 
of  the  melt.  The  crystal  is  allowed  to  rotate  at  an  angular  velocity  0,2  different  from  the  one 
of  the  crucible  Qi.  The  solid  vertical  and  bottom  boundaries  of  the  crucible  are  maintained 
^  at  constant  non-dimensional  temperature  of  Ti  =  2,  whereas  the  crystal-melt  interface  is  at 
a  fixed  temperature  of  72  =  1-  The  top  free  surface  is  assumed  to  have  negligible  heat  flux  to 
the  surrounding  gas  (an  assumption  which  is  sometimes  removed).  The  equations  of  motion 
are  the  incompressible  Navier-Stokes  equations,  together  with  the  Boussinesq  approximation 
for  the  effect  of  the  heating.  The  length  used  for  non-dimensionalization  is  the  radius  of  the 
crucible  R,  whereas  the  non-dimensionalizing  temperature  is  the  difference  Ti  —T2. 

The  boundary  condition  at  the  top  is  modified  in  order  to  include  radiation  or  other  kinds 
of  heat-loss,  whereas  the  base  of  the  crucible  can  be  insulated  instead  of  being  at  constant 
temperature.  The  process  of  crystal  growth  could  potentially  be  represented  by  a  suction 
velocity  Vc  at  the  crystal-melt  interface,  but  this  velocity  is  usually  negligible  when  compared 
to  the  motion  of  the  fluid. 

Assuming  that  we  always  work  in  the  frame  of  reference  of  the  rotating  crucible,  and 
that  the  Boussinesq  approximation  is  valid,  the  equations  of  motion  are 


(1) 


Depending  on  the  choice  of  the  reference  velocity  scale  U,  the  relevant  non-dimensional 
parameters  can  be  quite  different. 

For  a  problem  without  rotation  of  any  sort,  governed  purely  by  natural  convection,  and 
cissuming  that  the  geometry  is  fixed  and  the  aspect  ratios  H/D  (where  D  =  2R)  and  r/R 
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are  constant,  the  relevant  parameters  are  the  Grcisshof  number,  and  the  Prandtl  number, 
defined  as: 


Gr 


^gjTi  -  T2)R^ 


Pr 


V 
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and  the  corresponding  non-dimensionalized  equations  become: 


(2) 

(3) 


dv 

^  +  v.Vv  = 
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Kc 
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(4) 


when  U  is  chosen  to  be  (7  =  vfRGr^l'^.  Here  Re  =  Gr^!'^.  These  are  the  actual  forms  of  the 
equations  that  we  use  for  the  case  of  pure  natural  convection. 

In  the  case  of  rotation,  the  non-dimensionalization  can  be  based  on  either  a  diffusive  velocity 
scale  (as  in  the  non- rotating  case),  or  a  scale  which  is  defined  by  rotation,  i.e.  U  =  fli?. 
In  this  case,  there  is  one  more  non-dimensional  parameter  in  addition  to  Gr  and  Pr,  the 
Ekman  number  E,  defined  as: 


Assuming  that  U  is  still  U  =  v 
become: 


u 


(5) 


/RGr^^^,  the  corresponding  non-dimensionalized  equations 
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where  k  is  the  direction  of  rotation  0i  (here  aligned  with  the  z  axis),  and  Re  is  again  equal 
to  Re  = 

Finally,  when  there  is  also  differential  rotation  between  the  crucible  and  the  crystal  $72, 
the  Rossby  number  is  also  important: 


_  2(0i  -  O2) 
(0i  +  $12) 


(7) 


and  its  influence  comes  only  through  the  boundary  conditions  at  the  crystal-melt  interface. 


Following  Jones  [2]  the  range  of  parameters  may  be  subdivided  into  the  following  five 
regimes. 


Regime  I:  Heating. 

No  rotation  is  assumed.  The  relevant  problem  is  convection  in  a  vertically  heated  slot. 
The  flow  is  typically  three-dimensional  and  unsteady,  because  of  natural  convection,  for  most 
practical  ranges  of  Gr  number.  In  addition,  at  the  free  surface,  “spoke”  instability  patterns 
could  be  formed  due  to  heat  loss  to  surrounding  environment. 

Regime  II:  Crucible  rotation. 

In  pure  form  {Gr  =  0)  the  relevant  problems  are  stability  of  the  flow  around  a  rotating 
disk  (Ekman  layer)  and  the  stability  of  the  flow  inside  the  rotating  crucible,  [3],  [4]. 

Regime  III:  Differential  rotation. 

This  regime  includes  instabilities  connected  with  regime  II,  but  also  includes  new  insta¬ 
bilities  due  to  the  dynamics  of  the  shear  layer  between  the  rotating  crystal  and  the  crucible 
melt.  An  interesting  problem  to  study  is  at  which  Rossby  number  the  flow  becomes  non 
axisymmetric. 


Regime  IV:  Disk  rotation  with  heating. 

This  is  a  regime  where  geostrophic  motion  can  be  expected.  Usually  rotation  somewhat 
suppresses  instabilities,  even  though  unsteady  three-dimensional  flow,  because  of  baroclinic 
instabilities,  can  be  expected.  An  interesting  problem  to  study  is  what  is  the  optimal  rotation 


6 


frequency  that  maximally  suppresses  thermal  convection  instabilities. 

Regime  V:  Differential  rotation  with  heating. 

This  regime  is  similar  to  regime  IV,  but  the  motion  can  be  non  axisymmetric  because  of 
both  baroclinic  instabilities,  and  because  of  shear  layer  instabilities  arising  from  the  differ¬ 
ential  rotation. 


3  Numerical  Results:  Simulation  of  Melt  Flows 

The  numerical  experiments  that  we  have  performed  are  both  two-dimensional  and  three- 
dimensional  and  extend  all  the  way  from  Regime  I  with  only  natural  convection,  to  Regime 
V  which  includes  heating  and  differential  rotation.  Here  we  confine  ourselves  to  some  rep¬ 
resentative  results,  with  detailed  results  presented  in  papers  now  in  the  publication  process, 
[1].  The  code  used  for  all  these  studies  is  a  general  geometry  spectral  element  code  with 
high  order  accuracy  in  space  and  time.  Numerous  resolution  studies  (discussed  in  the  papers 
mentioned  above)  have  been  performed  to  ensure  the  accuracy  and  reliability  of  the  results 
presented  here. 

3.1  Crystal  melt  flow  -  2D  simulations 

3.1.1  High  aspect  ratio  crucible  H/D=l 

The  values  of  the  parameters  used  in  these  simulations  were  chosen  so  that  they  correspond 
to  the  values  used  by  [5]  in  the  simulation  of  buoyancy  driven  convection  in  Si  melts.  The 
particular  case  picked  was  one  in  the  unsteady  regime,  as  reported  in  [5],  and  corresponds 
to  a  value  of  Gr  =  2.8  x  10®,  with  Pr  —  0.03.  The  simulations  reported  in  that  paper  are 
three-dimensional  however.  The  aspect  ratio  of  the  crucible  for  this  set  of  simulations  is 
H/D  =  1. 

The  flow  starts  from  rest  and  Gr  is  gradually  increased  from  10, 000  to  40, 000  and  then  to 
2.8  X  10®.  The  first  two  simulations  reach  a  steady  state  and  then  the  flow  and  temperature 
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field  at  Gr  =  40,000  is  used  as  an  initial  condition  for  the  calculation  at  Gr  =  2.8  x  10®. 
Time  history  of  several  flow  variables  for  this  run  is  plotted  in  Figure  (2).  It  can  clearly 
be  seen  that  even  for  this  axisymmetric  case,  the  flow  is  unsteady  with  moderate  amplitude 
fluctuations.  Taking  into  account  that  the  axisymmetric  mode  is  usually  very  stable,  and 
that  it  is  typically  the  first  or  second  azimuthal  mode  that  goes  unstable  first,  this  result 
can  be  interpreted  as  an  instability  of  the  three-dimensional  flow  as  well. 

After  the  pure  boyancy  driven  flow  reaches  a  quasi  steady  state,  rotation  is  introduced 
in  order  to  study  its  effect  on  the  stabilization  of  the  flow.  The  magnitude  of  rotation 
introduced  gives  a  value  of  the  Ekman  number  E  =  Gr~^^^  ~  1.67  x  10“®,  which  is  typical  of 
rotation  frequencies  in  practice.  The  effect  of  rotation  on  the  flow  is  quite  significant,  since 
it  immediately  reduces  the  amplitude  of  the  fluctuations  as  seen  in  Figure  (2),  for  times  t  > 
87.5.  The  fact  that  the  flow  is  becoming  stabilized  with  rotation  in  the  axisymmetric  regime, 
does  not  necessarily  mean  that  it  is  also  stable  if  the  flow  is  allowed  to  be  three-dimensional. 
Therefore  three-dimensional  simulations  were  performed.  The  three  dimensional  simulations 
were  initiated  by  using  the  two  dimensional  field  with  rotation  at  t  =  172.5  as  an  initial 
condition. 

The  difference  between  the  flow  structure  with  and  without  rotation  is  illustrated  in 
Figures  (4),  and  (5),  respectively.  As  can  be  observed  in  the  Figures,  the  flow  with  no 
rotation  consists  of  the  recirculating  cells,  and  the  cold  fluid  is  able  to  convect  all  the  way 
to  the  bottom  of  the  crucible;  on  the  other  hand,  in  the  case  of  rotation,  there  exists  only 
one  recirculation  cell  and  the  amplitude  of  the  fluctuations  is  much  smaller.  The  difference 
between  the  amplitude  of  the  vorticity  field  generated  by  the  thermal  convection  for  the  two 
cases  is  demonstrated  in  Figure  (6)  where  the  maximum  vorticity  in  the  bulk  of  the  flow  is  of 
order  0(10)  without  rotation  and  0(1)  with  rotation.  For  both  cases  the  global  maxima  of 
vorticity  occur  on  the  side  walls  of  the  crucible,  as  expected  for  such  high  Grasshof  numbers. 

3.1.2  Low  aspect  ratio  crucible  H/D=0.25 

Several  simulations  were  carried  out  for  a  crucible  with  an  aspect  ratio  of  H/D  =  0.25, 
since  low  aspect  ratio  crucibles  are  more  common  in  practice.  The  runs  performed  for  this 
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Figure  4:  Instantaneous  flow  fleld  without  rotation  for  crucible  with  H/D=l,  axisymmetric 
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Figure  5:  Instantaneous  flow  field  with  rotation  for  crucible  with  H/D=l,  axisymmetric 
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Figure  6:  Instantaneous  vorticity  field  (a)  without  rotation  and  (b)  with  rotation,  for  crucible 
with  H/D=l,  axisymmetric 
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configuration  again  range  from  2-D  simulations  to  3-D  simulations  with  rotation  and  natural 
convection;  the  effect  of  crystal  differential  rotation  has  not  been  investigated  yet.  The  value 
of  the  Grasshof  number  used  in  the  simulations  was  picked  in  a  way  such  that  a  comparison 
between  the  high  {HID  =  1)  and  the  low  {HfD  =  0.25)  aspect  ratio  crucibles  is  meaningful. 
Keeping  all  other  parameters  the  same,  the  length  scale  used  for  non-dimensionalization  can 
be  either  L  =  {HRY^^  OT  L  =  {HR?Y^^->  which  results  in  an  equivalent  Gr  number,  for  the 
low  aspect  ratio  crucible,  equal  to  2.8  x  10®  or  4.4  x  10®  respectively.  This  value  is  equal 
or  higher  than  the  one  used  for  the  high  aspect  ratio  crucible  and  one  might  expect  a  more 
chaotic  flow  to  be  present.  In  addition,  one  has  to  take  into  account  the  fact  that  now  the 
total  volume  of  the  liquid  metal  is  double  the  one  corresponding  to  the  high  aspect  ratio 
crucible,  for  the  same  diameter  crystal  grown. 

The  2-D  simulations  at  the  Gr  number  mentioned  in  the  previous  paragraph,  reached  a 
steady  state  for  the  low  aspect  ratio  crucible.  Since  the  2-D  flow  for  the  high  aspect  ratio 
crucible  HfD  =  1  is  time  dependent  with  an  order  1  amplitude  of  velocity  fluctuations,  this 
fact  already  implies  that  flow  in  the  low  aspect  ratio  crucible  is  less  unstable  and  will  possibly 
lead  to  smaller  fluctuations  in  the  3-D  flow  as  well.  Moreover,  isocontours  of  azimuthal 
vorticity  indicate  the  presence  of  only  one  vortical  cell  structure  in  the  crucible,  whereas 
two  such  structures  are  present  in  the  axisymmetric  simulations  for  the  high  aspect  ratio 
crucible.  This  can  be  observed  in  Figure  (8)  where  isocontours  of  flow  variables  from  the 
steady  state  field  are  plotted.  As  can  be  observed  in  the  figure,  there  exists  a  small  pocket  of 
stagnant  fluid  close  to  the  lower  right  hand  corner  of  the  crucible,  which  could  be  a  potential 
source  of  instability;  this  implies  that  a  rounded  crucible,  following  the  streamlines  could  be 
more  effective  in  stabilizing  the  flow.  Time  histories  of  flow  components  at  several  points 
inside  the  crucible  are  plotted  in  figure  (7),  and  as  can  be  observed  all  values  reach  a  steady 
state  for  long  times. 

3.1.3  Rounded  low  aspect  ratio  crucible  H/D=0.25 

Besides  a  low  aspect  ratio,  the  crucibles  used  in  practice  are  rounded  at  the  bottom  so  that 
they  have  no  abrupt  corners.  Therefore,  our  current  and  near  future  studies  are  focusing  on 
crucible  shapes  similar  to  the  one  shown  in  figure  (9)  which  is  a  combination  of  a  cylindrical 
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part  on  top  and  an  ellipsoid  part  at  the  bottom.  The  simulations  performed  up  to  now  are 
only  axisymmetric  and  aa  can  be  observed  in  figure  (10),  the  flow  at  Gr  =  2.8  x  10®  reaches 
a  steady  state,  in  a  way  similar  to  the  flow  in  the  low  aspect  ratio  crucible  with  corners.  As 
can  be  observed  in  Figure  (9),  the  flow  again  consists  of  only  one  recirculating  cell,  but  now 
the  crucible  rounded  wall  follows  closely  the  flow  streamline  structure. 

3.2  Crystal  melt  flow  -  3D  simulations 

These  simulations  start  with  initial  conditions  from  the  2D  simulations  at  Gr  —  2.8  x  10®  by 
imposing  a  perturbation  on  the  first  azimuthal  mode  of  total  energy  of  the  order  of  10“^.  Two 
different  types  of  simulations  are  illustrated  here,  with  high  aspect  ratio  crucible  Hj D  =  1, 
and  with  low  aspect  ratio  Hj D  =  0.25. 

3.2.1  High  aspect  ratio  crucible  H/D=l 

1.  Heating,  Crucible  Rotation,  Differential  Crystal  Rotation 

An  example  is  simulation  with  rotation  of  both  the  crucible  and  differential  rotation  of 
the  crystal  with  Rossby  number  Ro  =  2.  These  simulations  show  that  the  flow  gradually 
develops  three-dimensionality  in  the  m  =  2  azimuthal  mode,  possibly  due  to  shear  layer 
instabilities  because  of  the  differential  rotation,  whereas  the  m  =  1  perturbation  decays  in 
time,  implying  that  baroclinic  instability  might  not  be  present.  In  general,  the  transition  to 
three-dimensionality,  is  a  much  slower  one  (than  the  natural  convection  case),  and  it  involves 
only  even  azimuthal  modes  (m=2,4,6,..).  The  tiihe  history  of  the  energy  of  the  first  three 
even  azimuthal  modes  is  shown  with  a  solid  line  in  Figure  (11);  the  first  three  odd  modes 
are  plotted  as  solid  lines  in  Figure  (15).  A  typical  azimuthal  energy  spectrum  is  shown  in 
Figure  (16)  where  the  difference  between  the  amplitudes  of  the  even  versus  odd  modes  is 
evident. 

Isocontours  of  instantaneous  flow  variables  on  both  r  —  z  and  r  —  (f)  planes  are  shown  in 
Figures  (12),  and  (13),  respectively,  and  eis  can  be  observed  -  especially  from  temperature 
isocontours  -  the  flow  demonstrates  an  even  m  =  2  symmetry. 
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Axisymmetric  Steady  state 

Gr=2.8xl0‘^ 


Axial  U:  -0.78, 0.96  Radial  V:  -1.01, 1.02 


Pressure  P:  -1.42, 0.065  Temperature  T:  0, 1 


Vorticity:  -8.6, 63.5 


Figure  9:  Isocontours  of  velocity,  pressure,  temperature,  and  vorticity  shown  for  rounded 
crucible  with  heating  at  Gr  =  2.8  x  10®. 
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The  structure  of  the  flow  in  terms  of  axial  vorticity  isocontours  is  shown  in  Figure  (14); 
as  can  be  observed  from  the  figure,  there  is  clearly  a  symmetry  in  the  flow  associated  with 
the  m  =  2  mode,  which  is  dominant,  demonstrating  itself  as  one  main  cyclonic  vortex  and 
two  counter-rotating  ones.  This  picture  is  similar  to  experimental  observations  reported  in 
[8],  where  the  experiment  involved  isothermal  rotation  of  a  disk  inside  a  rotating  crucible. 

2.  Heating 

Simulations  with  pure  natural  convection.  At  t  =  80.0,  and  after  the  three-dimensional 
flow  discussed  above  developed  some  three-dimensionality  showing  that  only  even  modes  are 
growing,  rotation  was  stopped  in  order  to  see  what  the  corresponding  final  states  of  the  flow 
would  be  without  rotation.  The  results  are  shown  as  dotted  lines  which  start  at  t  =  80.0,  in 
Figures  (11)  and  (15).  As  soon  as  rotation  is  stopped,  all  non-axisymmetric  modes  increase 
sharply  in  energy.  This  might  be  an  indication  that  for  the  range  of  parameters  used  in 
practice,  pure  natural  convection  can  be  violently  three-dimensional  and  unsteady,  whereas 
rotation  stabilizes  the  flow  and  reduces  the  fluctuation  levels.  In  some  cases  it  might  even  be 
that  rotation  prevents  three-dimensionality.  As  can  be  observed  in  Figures  (17),  and  (18), 
which  show  isocontours  of  flow  variables  on  r  —  and  r  —  ^  planes,  the  flow  is  quite  chaotic 
and  no  particular  symmetries  are  evident.  Also,  although  the  Prandtl  number  is  very  low, 
low  temperature  fluid  from  the  top  of  the  crucible  convects  all  the  way  to  the  bottom  of  the 
crucible,  indicative  of  the  large  amplitude  fluctuations  present  in  the  flow. 

3.  Heating,  Crucible  Rotation 

The  results  of  simulations  with  natural  convection  and  only  crucible  rotation  show  that  as 
soon  as  differential  rotation  of  the  crystal  is  stopped,  at  t  =  232.5,  the  three-dimensionality 
rapidly  increases,  since  now  all  non-axisymmetric  modes  (even  and  odd)  grow  in  time.  This 
can  again  be  observed  in  Figures  (11)-  (15)  by  noting  the  dotted  lines  which  start  att  =  232.5 
and  increase  in  amplitude.  The  amplitude  of  observed  fluctuations  in  this  case  seems  to  be 
much  larger  than  for  the  case  of  differential  rotation,  but  it  still  seems  to  be  smaller  than  for 
the  pure  natural  convection  case.  In  addition,  the  flow  structure  for  the  case  with  rotation 
(but  no  differential  rotation),  shown  in  Figures  (19)  and  (20),  is  not  as  chaotic  as  the 
pure  natural  convection  case,  described  in  the  previous  paragraph,  although  there  is  still  no 
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Azimuthal  even  mode  time  history 


Figure  11:  Time  history  of  even  azimuthal  mode  energies 


Rotation— Diff.  Rotation 

Gr=2.8xl0\  E=6xl0^,  Ro=l 


Azimuthal  W  Temperature  T 


Figure  12:  Instantaneous  isocontours  of  flow  variables  on  a  r  —  2:  plane  for  case  with  heating 
rotation,  and  differential  rotation 
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Rotation— Diff.  Rotation 

Gr=2.8xl0\  E=6xlO~tRo=l 


Z=1 


Azimuthal  W 


z=l 


Radial  V 
z=L75 


Temperature  T 


Figure  13:  Instantaneous  isocontours  of  flow  variables  on  a  r  —  ^  plane  for  case  with  heating 
rotation,  and  differential  rotation 


23 


Figure  14;  Isocontours  of  axial  vorticity  at  z  =  1  for  case  with  heating  rotation,  and  differ¬ 
ential  rotation 
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No  Rotation 

Gr=2.8xl(f 


Azimuthal  W 


Temperature  T 


Figure  17:  Instantaneous  isocontours  of  flow  variables  on  a  r  —  2:  plane  for  case  with  only 
heating 
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No  Rotation 

Gr=2.8xl(f 


Z=1  z=l 


Axial  U  Radial  V 


z=l 


Azimuthal  W 


z=1.75 


Temperature  T 


Figure  18:  Instantaneous  isocontours  of  flow  variables  on  a  r  —  ^  plane  for  case  with  only 
heating 
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evidence  of  any  synunetries  present  in  the  flow. 


3.2.2  Low  aspect  ratio  crucible  H/D=0.25 

The  3-D  runs  were  performed  using  the  axisymmetric  steady  state  result  as  initial  condition, 
with  a  10“®  perturbation  in  the  m  =  1  azimuthal  mode.  As  figure  (21)  shows,  the  even  modes 
(m  =  2,4...)  increase  exponentially  in  amplitude,  whereas  the  odd  modes  don’t  start  to  pick 
up  until  after  the  even  modes  have  saturated  to  a  certain  amplitude.  Rotation  at  an  angular 
velocity  Qj  =  1  is  turned  on  at  t  =  247.5.  The  flow  becomes  three  dimensional  for  both 
cases  (with  or  without  rotation),  however,  the  odd  modes  grow  faster  for  the  case  without 
rotation  and  all  modes  saturate  at  slightly  lower  amplitudes  when  rotation  is  turned  on.  A 
typical  instantaneous  flow  field,  shown  in  terms  of  velocity  and  temperature  isocontours,  is 
plotted  in  figures  (22),  and  (23),  for  the  case  with  both  heating  and  rotation.  It  can  be 
observed  in  these  figures  that  there  exists  a  dominance  of  even  modes  in  the  flow,  and  in 
particular  of  the  m  =  4  azimuthal  mode. 

3.2.3  Rounded  low  aspect  ratio  crucible  H/D=0.25 

Preliminary  3-D  simulations,  without  rotation,  indicate  that  the  flow  also  becomes  three- 
dimensional,  starting  with  the  even  modes,  but  these  simulations  have  to  be  continued 
further.  It  is  interesting  to  note  that  since  the  flow  in  low  aspect  ratio  rounded  crucibles 
seems  to  be  less  unstable  than  high  aspect  ratio  crucibles  with  corners,  one  might  be  able  to 
perform  full  three-dimensional  numerical  simulations  at  parameters  {Gr,  E,  Ro)  which  are 
close  to  realistic  values  used  in  practice;  that  is  to  say  that,  although  the  same  flow  without 
rotation  or  even  without  differential  rotation  might  be  too  chaotic  to  resolve  with  a  DNS, 
the  flow  with  all  the  stabilizing  mechanisms  turned  on  could  be  not  out  of  reach.  This  does 
not  necessarily  mean  that  the  resulting  flow  is  axisymmetric. 
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Rotation 

1/2 

Gr=2.8xl0f  E=Gr~  =6x10 


Azimuthal  W  Temperature  T 


Figure  19:  Instantaneous  isocontours  of  flow  variables  on  a  r  —  z  plane  for  case  with  heating 
and  rotation 
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Rotation 

1/2 

Gr=2.8xl0,  E=Gr~  =6x10 


z=l  z=l 


z=l 


Azimuthal  W 


z=L75 


Temperature  T 


Figure  20;  Instantaneous  isocontours  of  flow  variables  on  a  j —  (j)  plane  for  case  with  heating 
and  rotation 
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Rotation:  Low  aspect  Ratio,  Gr=2.8xlO^,  E=6xl0  ^ 


Figure  21:  Time  history  of  3-D  mode  energies  for  flow  in  low  aspect  ratio  crucible  at  Gr 
2.8  X  10® 
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I^w  aspect  ratio  crucible: 

Rotation-thermal  convection 

Gr=2,8xl0\  E=Gr~^'"=6xlO~^ 

z=0.5  z~0.5 


Axial  U 


Radial  V 


z=0.5 


Azimuthal  W 


z=0.875 


Temperature  T 


Figure  22:  Isocontours  of  axial  velocity  and  temperature  shown  for  a  top  view  for  the  low 
aspect  ratio  crucible  with  heating  and  rotation. 
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Tmw  aspect  ratio  crucible: 

Rotation— thermal  convection 

Gr=2.8xl0‘,  E=Gr~‘'^=6xlO^ 


Axial  U 


Azimuthal  W 


Temperature  T 


Figure  23;  Isocontours  of  axial  velocity  and  temperature  shown  for  a  front  view  for  the  low 
aspect  ratio  crucible  with  heating  and  rotation. 
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3.3  Magnetic  field  simulations 

Two  different  formulations  have  been  developed  and  implemented  for  the  case  of  an  external 
magnetic  field  applied  to  the  crystal  melt.  In  the  first  one,  we  solve  the  full  magnetic  field 
and  hydrodynamic  equations,  whereas  in  the  other  we  assume  that  the  magnetic  Reynolds 
number  Rcm  is  very  small  thereby  allowing  simplification  of  the  solution  procedure. 

The  governing  MHD  equations  at  the  low  magnetic  Reynolds  number  limit  are 

—  +  v  •  Vv  =  -Vp - — —  +  ^  W  +  r  +  TVj  X  Bo 

at  Re  E  Re 

j  =  -V^  + V  X  Bo 


VV  =  V  •  (v  X  Bo) 

|^  =  n.(vxB„) 

Here,  Re  is  again  defined  as  Re  =  Bo  is  the  imposed  axial  magnetic  field,  and  <f>  is 

the  electric  potential  corresponding  to  the  induced  electric  field.  We  present  our  results  in 
terms  of  the  magnetic  parameter  N  =  -^Rem-  N  is  related  to  the  Hartmann  number  by 

R2 

Ha'^  =  NRe  =  —^Re^Re 
pU^ 

Various  simulations  have  been  performed  with  the  main  focus  being  to  understand  the  in¬ 
fluence  of  the  magnetic  field  on  the  melt  flow.  All  simulations  performed  up  to  now  employ 
an  axial  magnetic  field,  since  a  transverse  one  would  automatically  destroy  axisymmetry 
(although  will  not  necessarily  increase  fluctuations).  From  our  simulations,  it  seems  that  the 
imposed  magnetic  field,  depending  on  its  amplitude,  can  either  dampen  three  dimensionality 
or  enhance  it.  This  can  be  observed  in  Figure  (24),  where  one  sees  that  a  magnetic  field 
with  a  magnetic  parameter  N  =  \  dumps  all  non-axisymmetric  modes,  whereas  a  value  of 
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20  increases  the  amplitude  of  non-axisymmetric  modes.  In  fact,  an  intermediate  value,  e.g. 
N  =  b,  reduces  the  amplitude  of  the  non-axisymmetric  modes  but  not  to  zero.  In  order  to 
visualize  the  effect  of  the  magnetic  field  on  the  flow  field,  stream-wise  vorticity  isocontours 
at  z  =  1  are  plotted  in  figures  (25a),  (25b),  (25c)  and  (25d).  As  can  be  seen  in  these 
figures,  the  resulting  vortex  shape  for  A  =  1  (Figure  25a)  is  a  circle,  corresponding  to  an 
axisymmetric  flow;  on  the  other  hand,  for  A  =  5  (Figure  25b),  the  vortex  shape  is  only 
slightly  distorted  from  axisymmetry,  whereas  for  A  =  10  (Figure  25c)  and  20  (Figure  25d) 
^  the  main  central  vortex  is  distorted  from  a  circular  shape  to  an  ellipse  with  a  higher  aspect 
ratio  than  the  case  without  magnetic  field,  shown  in  figure  (14).  In  addition,  the  vortex 
structure  rotates  with  respect  to  the  rotating  crucible. 

We  have  also  extended  the  analysis  due  originally  to  Chandrasekhar  (Oxford  University 
Press,  1961)  to  analyze  the  effects  of  magnetic  fields  on  the  stability  of  Benard  convection. 
This  analysis  which  involves  the  solution  of  a  12th  order  boundary  value  problem,  leads 
to  interesting  insights  on  the  mechanisms  by  which  magnetic  fields  can  be  destabilizing,  as 
shown  in  the  following  section. 

3.3,1  Influence  of  rotation  and  magnetic  field  on  the  stability  of  thermal  con¬ 
vection 

The  growing  crystall  is  usually  rotated  as  it  is  pulled.  The  objective  is  to  improve  uniformity 
by  providing  a  viscous  shear  layer  that  tends  to  isolate  the  growth  interface  from  the  turmoil 
deeper  in  the  melt.  The  crucible  is  also  rotated  to  smooth  out  thermal  asymmetries  that 
might  arise  from  irregularities  in  the  heating.  Sometimes  the  crystal  melt  is  stabilized  by 
the  use  of  a  magnetic  field.  Crystal  rotation  and  magnetic  field  individually  may  suppress 
curtain  instabilities,  but  combined  together  they  may  lead  to  new  instabilities. 

As  an  example  of  new  instabilities,  which  could  potentially  arise  by  the  combined  use 
of  rotation  and  magnetic  field,  we  consider,  following  Chandrasekhar  (1961)  [9],  the  case  of 
Rayleigh-Benard  convection  in  the  presence  of  rotation  and  magnetic  field.  The  convection 
is  characterized  by  the  Rayleigh  number  Ra,  Prandtl  number  Pr,  Taylor  T  or  Ekman  E 
numbers,  with  the  Taylor  number  being  T  =  1/A^,  and  the  Hartmann  number  Ha. 
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MHD:  Even  azimuthal  modal  energy  in  time 


Figure  24:  Time  history  of  the  even  azimuthal  mode  energy  with  the  application  of  different 
magnetic  fields 
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c)  N=10  d)  N=20 


Figure  25:  Isocontours  of  axial  vorticity  at  2  =  1  for  the  case  with  a  magnetic  field  with  a) 
iV  =  1,  b)  N=5,  c)  N=10,  and  d)  N=20 
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Separately  both  rotation  and  magnetic  field  inhibit  the  onset  of  instability  and  they  both 
elongate  the  cells  which  appear  at  marginal  stability.  These  effects  have  a  common  origin: 
the  flow  becomes  more  two-dimensional.  Acting  together,  however,  they  may  lead  to  a  new 
instability.  The  reason  behind  this  is  that  rotation  and  magnetic  field  present  together  have 
conflicting  interests.  Rotation  induces  a  component  of  vorticity  in  the  direction  of  rotation. 
It  results  in  the  streamlines  becoming  closely  wound  spirals  with  motions  principally  confined 
to  planes  transverse  to  the  rotation.  On  the  other  hand,  a  magnetic  field  would  suppress  the 
motion  transverse  to  its  direction  and  the  motion  along  the  magnetic  field  becomes  dominant. 
In  addition,  the  instability  in  the  presence  of  rotation  sets  in  as  overstability,  although  it 
sets  in  as  stationary  convection  when  a  magnetic  field  is  present. 

In  the  simplest  case  of  two  free  boundaries  the  characteristic  equation  for  the  critical 
Rayleigh  number  takes  the  form 

^  =  ^  - xl(l  +  xf  +  Ha?] 

where 


fi  rr 

x  =  and  Ti  =  —  (9) 

tt’  it 

and  i  is  the  characteristic  scale  of  the  onset  of  instability.  These  equations  determine  the 
instability  threshold  in  the  case  of  the  onset  of  instability  as  stationary  convection.  In  the 
case  of  the  onset  of  instability  as  overstability  the  characteristic  equations  take  the  form 


X 


(10) 


where 


2  _  (1  +  -  Pr)  -  PrHa'^ 

^  1  +  a:  (1  -|-  a;)2(l  -|-  Pr)  +  PrHa? 

These  equations  are  accurate  in  the  case  of  low  magnetic  Prandtl  number  which  is  always 
the  case  for  crystal  growth. 


1  +  a:  -t- 


Ha^ 


1  -f  X 


(11) 
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The  solution  of  these  equations  for  Ekman  number  equal  to  =  10  ^  (or  Taylor  number 
10®)  and  Prandtl  number  Pr  =  0.2  is  shown  in  figure  (26).  Let  us  consider  the  stability  of 
the  flow  with  the  increase  of  the  magnetic  field  starting  from  Rayleigh  numbers  for  which 
the  flow  is  stable  in  the  absence  of  the  magnetic  field.  As  can  be  observed  the  flow  may 
become  unstable  at  a  certain  level  of  the  magnetic  field.  If  we  increase  the  strength  of  the 
magnetic  field  even  further  the  flow  stabilizes  once  again. 

The  effect  of  the  destabilization  of  thermal  convection  by  the  action  of  a  magnetic  field  in 
the  presence  of  rotation  is  sensitive  to  the  geometry,  Prandtl  number  and  configuration  of 
temperature  gradients.  The  stability  of  the  flow  in  crystal  melts  in  the  presence  of  magnetic 
field  and  rotation  should  be  studied  using  three-dimensional  direct  numerical  simulation. 


4  RNG-Based  Modeling  of  Stratified  Mixing 


Renormalization  group  (RNG)  methods  are  a  general  framework  for  “model  building”  in 
which  the  complex  dynamics  of  physical  problems  is  described  in  terms  of  so-called  “coarse¬ 
grained”  equations  of  motion  governing  the  large-scale,  long-time  behavior  of  the  physical 
system.  The  RNG  approach  allows  coarse- graining  of  physical  phenomena  as  varied  as  critical 
phenomena,  high-energy  physics,  and,  especially  in  the  context  of  fluid  dynamics,  turbulence, 
combustion,  and  heat  transfer.  The  key  idea  is  that  the  RNG  method  is  applicable  to 
scale  invariant  phenomena  lacking  externally  imposed  characteristic  length  and  time  scales. 
For  turbulence,  this  means  that  the  method  is  applicable  to  the  description  of  the  small 
scales  (small  eddies)  that  should  be  statistically  independent  of  the  external  conditions  and 
dynamical  forces  that  create  them  through  various  kinds  of  instability  phenomena.  In  other 
words,  the  RNG  method  gives  a  theory  of  the  so-called  Kolmogorov  equilibrium  range  of 
turbulence,  especially  comprising  the  so-called  inertial  range  of  small-scale  eddies  whose 
energy  spectrum  follows  the  famous  Kolmogorov  law  E{k)  ~  The  importance  of  the 

RNG  results  is,  as  we  shall  see,  once  the  inertial  range  eddies  can  be  accounted  for  in  a 
qualitatively  correct  way,  we  may  then  obtain  coarse-grained  equations  of  motion  for  the 
other  relevant  variables  of  the  turbulence,  including,  the  mean  velocity,  rms  velocities,  etc. 
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stability  curve  for  E  =  10“^^,  Pr  =  0.2 


logic  Ha 

Figure  26:  The  stability  diagram  of  critical  Rayleigh  number  as  the  function  of  Hartmann 
number  for  Ekman  number  10~^ 

Summary  of  the  RNG  Method 

In  the  application  of  RNG  methods  for  turbulence  modeling,  the  local  turbulent  kinetic 
energy  K  and  the  local  energy  dissipation  rate  £  are  used  to  eliminate  the  large-eddy  length 
scale  L  from  the  dynamics,  and  the  RANS  equation  of  motion  for  Ui  is  supplemented  by 
equations  for  K  and  £.  The  RNG  method  is  used  to  evaluate  otherwise  unknown  terms  in 
tehse  equations  for  this  K  —  £  model. 

The  resulting  RNG  K  —  £  model  differs  from  classical  (or  standard)  K  —  £  models  in  at 
least  six  ways; 
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•  At  high  Reynolds  number,  the  constants  in  the  RNG  K  -  £  model  are  evaluated  by 
the  theory; 

•  New  terms  appear  like  a  rate  of  strain  term,  which  is  important  for  treatment  of  non¬ 
equilibrium  effects  and  flows  in  the  rapid  distortion  limit; 

•  Impulse  response  modifications,  important  for  non-equilibrium  effects,  are  also  taken 
into  account; 

•  Low  Reynolds  number  modifications  are  given  by  the  RNG  theory,  so  wall  function  are 
no  longer  required; 

•  Modified  boundary  conditions  are  developed  from  the  theory; 

•  Stratification  and  rotation  (swirl)  effects  are  accounted  for  by  extensions  of  the  RNG 
theory  to  stratified  and  rotating  flow.  Here  we  briefly  describe  these  new  effects  as 
developed  in  the  present  work. 

RNG  modeling  of  stratified  turbulence 

We  have  initiated  the  derivation  of  turbulence  transport  models  for  vertical  mixing  based 
on  the  RNG  analysis  of  stratified  flows.  This  model  will  be  used  for  the  investigation  of  highly 
turbulent  flow  inside  the  stratified  gas  surrounding  the  crystal  and  the  melt,  or  for  very  high 
Rayleigh  numbers  for  the  liquid  flow  inside  the  crucible  as  well.  The  classical  theory  of 
turbulent  thermal  convection  is  based  on  the  following  main  ideas: 

1.  The  physics  of  boundary  layers  is  independent  of  the  processes  taking  place  in  the  bulk 
of  the  fluid  outside  the  thermal  boundary  layers. 

2.  The  energy  and  temperature  spectra  are  E{k)  ~  ET{k)  ~  k~^l^.  These  assumptions  are 
sufficient  to  derive  the  classical  result  that  the  Nusselt  number  Nu^  a  non-dimensional 
heat  flux  ,  satisfies 


Nu  ~  Ra^'^ 
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where  Ra  is  the  Rayleigh  number.  In  recent  years,  these  classical  assumptions  have  been 
brought  into  question  by  Libchaber’s  experiments  at  Chicago  which  covered  a  broad  range 
of  Rayleigh  numbers  from  10’’  to  more  than  10’^.  The  surprising  result  of  these  carefully 
controlled  experiments  was  that  E-riJt)  ~  and  Nu  ~  . 

Recently,  a  theoretical  explanation  of  these  experiments  has  been  given  by  V.  Yakhot 
at  CHI.  The  theoretical  result  shows  that  at  intermediate  Rayleigh  numbers  the  Nusselt 
number  scales  as  while  at  very  large  Rayleigh  numbers  Nu  ~  At  Rayleigh 

numbers  of  order  of  10^^,  the  classical  Ra}^^  theory  predicts  heat  fluxes  that  are  an  order 
of  magnitude  larger  than  those  found  by  either  the  Libchaber  experiments  or  the  Yakhot 
theory. 

An  important  outcome  of  the  new  theory  is  that  the  large-scale  flow  dynamics  is  governed 
by  not  the  energy  flux  S,  but  rather  by  entropy  flux  N  at  high  Rayleigh  numbers.  Here,  N 
is  given  by  Y  =  K{VTy.  Consequently,  the  effective  viscosity  is  estimated  as: 


l/y 


c:. 


KT^ 

N 


(12) 


where  K  is  the  turbulence  kinetic  energy  and  (7^  is  a  numerical  constant.  This  formula 
replaces  the  K  —  S  model 


VT  =  - , 

e 

for  turbulent  flows  dominated  by  the  energy  flux.  Therefore,  in  order  to  successfully  model 
convective  turbulence,  equations  for  and  N  are  needed.  The  renormalization  group 

analysis  applied  to  this  problem  leads  to  a  set  of  equations  for  and  These 

K  —  £  —  N  turbulence  transport  equations  have  a  form  analogous  to  that  of  the  RNG  K  —  £ 
equations  which  have  already  been  successfully  applied  to  complex  turbulent  flows  with 
strong  rotation,  massive  separation,  flow  transitions  and  other  phenomena. 

RNG  description  of  anisotropic  transport  in  stratified  turbulence 

The  system  of  equations  for  mean  velocity,  f/,,  mean  temperature,  T  (relative  to  an 
’ambient’  temperature  To),  mean  pressure  T,  mean  turbulent  kinetic  energy,  K  and  mean 
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energy  dissipation  rate,  e  are 


dt 


'dxj 


m 

dxi 


=  0 


dK  „  ,vh^  9  r’^zdK,  , 

—  =  V;i( — WhK)  +  - — ( —  —  )+  „ 

dt  dxj  (Tk  0x3  akOxz  poOXj 


TijdUi  ,  dT 
- 9^1^. - 


de 


de 


^  +  c^iir  =  v.ev.e)  +  — ( 


d  ,Vz  de 


TijdUi 


dt 


dx 


dxs^a.dxz^'^  Po  dxj 


dxz 
dT_ 
3 


C pQI^^Z  r\ 

dxz 


Let  us  introduce  parameter  of  stratification  (buoayancy  frequency) 

dT, 


dx- 


a  =  /{gl3\ 

Corresponding  dimensionless  parameter  is 

C  =  nK/e 


■I) 


Following  Prandtl  numbers  can  be  defined  as 

QJl  =  ^h/f^isot 
OC2  =  J^hlVz, 
OCZ  =  Vhll^h, 
OL4  =  Vh/Hz 


Here  Uiso  =  All  Prandtl  numbers  are  functions  of  C,  only. 

The  system  of  the  RNG  equations  for  these  anisotropic  Prandtl  numbers  has  the  following 
form: 


=  /.(ai,  0:2, 0:3, 0:4,  C)- 


i  =  1,2, 3, 4 


where  the  RNG  theory  determines  the  explicit  form  of  the  function  /,.  It  is  these  equations 
that  provide  for  the  effective  treatment  of  stratified  flows. 
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5  Future  work 


The  work  described  here  provides  significant  insights  into  the  kinds  of  fluid  dynamical  effects 
important  to  the  MLEK  process.  In  particular,  our  results  on  flow  stabilization  by  rotation, 
differential  rotation  and  magnetic  fields  should  have  substantial  impact  on  the  choice  of 
operational  flow  regimes  for  successful  crystal  growth  by  the  MLEK  process.  However,  much 
additional  work  needs  to  be  done.  In  particular,  additional  work  on  the  following  kinds  of 
problems  is  important  for  practical  applications: 

•  Heat  losses  to  the  gas:  Detailed  consideration  of  modification  of  the  flux  boundary 
condition  at  the  free  surface  must  be  given  to  account  for  typical  radiative  or  convective 
heat  losses  to  the  surrounding  gases. 

•  Interface  tracking:  An  adaptive  remeshing  algorithm  for  interface  tracking  will  have 
significant  impact,  since  the  interface  shape  between  the  melt  and  the  crystal  and 
between  the  melt  and  the  gas  has  to  dynamically  be  updated  during  the  calculation. 
For  this,  a  spectral  element  adaptive  algorithm  similar  to  the  one  for  moving  boundaxies 
will,  we  believe,  be  most  effective. 

•  Thermal  stresses  in  crystal:  Having  calculated  the  shape  and  and  temperature  distri¬ 
bution  along  the  interface,  the  temperature  and  stress  distribution  inside  the  crystal 
can  be  obtained.  From  the  stress  distribution,  a  density  of  dislocations  can  be  obtained 
inside  the  crystal  which  can  give  information  about  the  quality  of  the  crystal  produced. 

•  Multi-species  composition:  Another  interesting  direction  is  the  incorporation  of  differ¬ 
ent  species  concentrations  as  part  of  the  melt  composition.  This  is  now  straightforward 
numerically,  but  its  difficulty  lies  in  the  fact  that  initial  and  boundary  conditions  might 
not  be  available  or  traceable  from  experiments. 
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